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Abstract 


V 

This  paper  describes  an  algorithm  and  a  FORTRAN  program  called 
MCHAIN  for  simulating  k  parallel  Monte  Carlo  replications  of  a  Markov 


V 

chain  using  rotation  sampling.  This  method  of  sampling  produces  k 
sample  transition  frequency  vectors  with  a  desirable  structure  of 
statistical  dependence  among  them.  In  particular,  these  sample  vectors 

I  yyi>  r 

can  be  used  to  estimate  the  probability  that,  say,  n^,...,ny  transitions 
of  types  l,...,r  occur  during  a  first  passage  from  state  a  to  state  b 

rh 

with  a  variance  of  the  estimate  of  0(1/K  f  and  a  computation  time  0(k) 
as  k  -+  »  .  This  compares  favorably  with  the  case  of  independent 

replications  wherein  the  estimate  would  have  a  variance  0(l/k)  and 

\ 

computation  time  0(k)  as  k  ->■  »  .  An  example  including  a  sample  driver 
program  are  presented  to  illustrate  how  MCHAIN  works  in  practice. 

.  I 
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Introduction 

Let  p  *  ||Pj|||  denote  a  positive  recurrent  aperiodic  Markov 
Chain  with  state  space  S  =  (0,l,...,n)  .  Let  N..  denote  the  number 
of  one-step  transitions  that  occur  from  state  i  to  state  j  during 

a  first  passage  from  state  a  to  state  b  a,b  e  S  on  an  arbitrarily 

(l) 

selected  replication  and  let  N. .v  '  denote  this  quantity  observed  on 

*  J 

replication  l  .  The  purpose  of  this  paper  is  to  describe  an  algorithm 
and  a  FORTRAN  program  called  MCHAIN  that  generates  k  sample  data  sets 
{N  0 ) .  i ,3  c  S},...,{N.  i,j  e  S}  from  the  chain  p  with  a  specialized 

dependence  among  sets  that  allows  one  to  estimate  quantities  of  interest 
with  greater  statistical  efficiency  than  independent  data  sets  allow. 

As  an  example  of  the  use  of  these  data  sets,  suppose  one  wants  to 


compute 

p(n^f...,n  )  =  probability  that  during  a  first  passage  (1) 

from  a  to  b  n.,...,n.  one-step  transitions 
occur  from  i^  to  j^»...,i  to  jr  ,  respectively. 
For  Markov  chains  with  relatively  general  structure,  convenient 
analytical  representations  are  not  available  for  (1)  so  that  computation 
is  not  possible.  In  fact  no  convenient  representations  are  available 
for  chains  with  special  structure  but  arbitrary  r  >  1  .  One  way  to 
overcome  this  absence  of  representation  is  to  employ  Monte  Carlo  methods. 
Assume  there  are  k  simulated  data  sets  {N..^;  i ,  j  e  S},..., 

*  sj 

Ik) 

{ N . . v  ' ;  i,j  e  S)  available  to  estimate  (1)  ,  each  simulation  beginning 

•  J 

with  a  departure  from  state  a  and  ending  upon  entry  into  state  b  .  Let 


Ci  5"v  j. 


(2) 
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where 


if  » 


u>. 


«(x)  =1  if  x  =  0 

=  0  otherwise  . 

i,j  £  S}  l  =  l,...,k  are  independent,  then 


(3) 


is  an  unbiased  estimator  of  p(n.| ,. . .  ,nr)  with  variance  proportional 

to  1/k  .  Using  a  specialized  sampling  procedure  called  rotation  sampling, 

ll) 

MCHAIN  induces  statistical  dependence  among  the  data  sets  {N..v  i,j  c  S) 

*  J 

Z  =  l,...,k  that  preserves  the  unbiasedness  of  (3)  but  changes  the  convergence 
rate  of  its  variance  to  0(l/k  )  .  Fishman  (1982b)  illustrates  the  significance 
of  this  improved  convergence  for  the  estimation  of  the  probability  density 
function  of  first  passage  time  in  a  semi-Markov  process.  More  generally, 
the  improved  convergence  is  clearly  beneficial  when  estimating  functions 
of  Markov  chains  that  are  linear  combinations  of  p(n^,...,nr)  for  varying 
n^,...,nr  .  Complete  details  about  rotation  sampling  can  be  found  in 
Fishman  (1982a). 

Section  1  describes  how  one  constructs  a  simulation  model  to  induce 
the  desired  dependence  and  shows  its  consequences  in  terms  of  the  improved 
efficiency.  For  comparison  a  procedure  is  first  described  for  simulating 
independent  replications  of  a  Markov  chain.  Section  2  describes  MCHAIN 
and,  in  particular,  its  input  requirements.  Section  3  presents  an  example 
of  how  MCHAIN  performs. 

Serial  and  Parallel  Simulation 

To  put  the  proposed  sampling  scheme  into  perspective  we  first  describe 
the  more  conventional  simulation  of  a  Markov  chain  in  which  data  on  k 
independent  replications  are  collected  in  serial  order. 


Serial  Simulation 

1.  Start  in  state  i=a  . 

2.  For  m=l . k  and  i,j=0,l . n:  N..^«-0. 

^  J 

3.  m  ■*-  1  . 

4.  Sample  U  from  U( 0,1)  . 

Use  U,  j=0,l,...,n}  and  the  alias  sampling  method 

(Walker  1977  and  Kronmal  and  Peterson  1979)  to  determine 
entered  state  j  . 

5  N  M  «.  N  +  1 

6.  If  j?b  ,  i  •*-  j  and  go  to  4  . 

7.  If  m=k,  deliver  (N  i,j=0,l,...,n}  m=l,...,k  . 

*  J 

8.  m  m+1 

9.  i  4-  a  . 

10.  Go  to  4  . 

Here  U(0,1)  in  step  4  denotes  the  uniform  distribution  on  the  half 
open  interval  [0,1)  .  The  alias  sampling  method  is  a  procedure  for 
sampling  from  a  tabled  discrete  distribution  at  a  cost  independent  of 
the  length  of  the  table.  For  example,  if  the  entered  state  is  iVb  , 
then  sampling  occurs  from  the  distribution  p.Q,  p^,...,p.  with 
constant  execution  time  regardless  of  the  value  of  n  . 

In  the  case  of  parallel  simulation  based  on  rotation  sampling  all 
k  replications  begin  at  the  same  moment  in  time.  The  steps  are: 
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Parallel  Simulation 

Definitions:  im  =  state  occupied  by  replication  m 

K.  =  number  of  replications  in  state  j  at  the 
3 

beginning  of  a  step. 

K*  =  number  of  replications  in  state  j  at  the 
J 

end  of  a  step. 


1.  For  m=l,...,k: 

la.  i  *■  a  . 

m 

lb.  For  i=0»l,...,n;  j=0,l,...,n  ,  0  . 

2.  For  i=0,l ,. . .  ,n  :  ■*-  0  . 

3.  For  i=0,l ,. ..  ,n  :  K*  ■*-  0  . 

4.  Sample  Ug,  U,,...,U  independently  from  U(0,1)  . 

5.  For  m=l , . . . ,k  : 

5a.  If  i  =b  ,  go  to  5h. 
m 

5b.  Determine  j  using  U.  and  (p.  .  ;  j=0,l,...,n} 

'm  mJ 

with  the  Fishinan-Moore  (1982)  cutpoint  sampling  method. 

5c.  K*  «-  K*  +  1  . 

J  J 


5d.  N.  , 
mJ 


(m) 


N.  +  1  . 

V 


5e. 


+  1/K, 


m 


m 


m 


5f .  If  U.  ^  1  ,  U.  <-  U.  -  1 


m 


m 


m 


5g.  im  -  j 


5h.  Continue  . 


K,  4-  Kk  +  Kf  . 
b  b  b 

If  K  =k  ,  deliver  {N. . 
b  i  j 


(m) 


^ »J”H|1 » •  * .  >n}  . 


8.  For  j=l,...,k  and  j^b  :  K.  *■ 


K* 

J 


9.  Go  to  3. 
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Here  rotation  sampling  is  implemented  in  step  5e  so  that  for  K. 

m 

replications  in  state  i  the  K.  uniform  deviates  used  to  determine 

m  i 

m 

the  newly  entered  states  have  the  uniform  distribution  on  [0,1)  but 
are  not  independent.  Fishman  (1981,  1982b)  contain  more  detailed 
descriptions  of  rotation  sampling.  Step  5b  uses  an  alternative  method 
to  the  alias  to  determine  the  entered  state.  In  particular,  this 
cutpoint  method  uses  the  inverse  transform  method  to  select  the 
entered  state  for  replication  m  as 


j  =  (r:  Pi  s  *  U.  <  j£=0  p.  s,  pi  ,-0  ;  r=0,l . n)  (4) 

mm  mm 

at  a  cost  independent  of  n  .  This  procedure  for  sampling  from  the 

Markov  chain  together  with  the  rotation  sampling  induces  the  desired 

dependence.  Since  the  alias  sampling  method,  which  is  slightly  less 

costly,  does  not  use  the  inverse  transform  method,  it  cannot  be  used 

to  induce  the  desired  correlation. 

To  appreciate  the  significance  of  the  dependence  that  rotation 

sampling  and  the  cutpoint  method  together  induce,  we  first  focus  on  (2). 

Let  w=n.+. . .+n^  .  Then 
1  r 


Xw"ni  ...  j;w"nT"'"nr-i  K 
n2=0  nr=0 


n 


1 


.n 


(5) 


is  the  number  of  replications  absorbed  on  the  wth  transition.  Making 
use  of  the  fact  that  var  Kw  =0(1)  ,  Proposition  1  in  Fishman  (1982b) 
shows  that  var  Kn^  n  =0(1)  .  Therefore,  var  p(n1,...,nr)  =  0(l/k2)  . 
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Computation  Time  Complexity 

Although  the  accelerated  convergence  of  the  variances  of  estimators 

is  important,  one  also  needs  to  consider  how  the  computation  time 

complexity  of  the  parallel  simulation  (PS)  compares  with  that  of  the 

serial  simulation  (SS).  Both  are  0(K)  .  Although  for  a  given  k  parallel 

simulation  takes  more  time  than  serial  simulation  does,  the  ratio  of 

these  times  converges  to  a  constant  as  k  ®  ,  thereby  demonstrating 

the  superiority  of  parallel  simulation  as  measured  by  the  common  variance 

reduction  measure  (see  Hammersley  and  Handscomb  1964) 

ud  -  variance  using  SS  .  computation  time  using  SS  _  n,.i  k  /fi\ 
variance  using  PS  computation  time  using  PS  '  '  ' 

2.  The  MCHAIN  Program 

Figure  1  lists  a  FORTRAN  subroutine  called  MCHAIN  that  generates 
the  transition  frequency  data  for  a  finite-state  Markov  chain  using 
rotation  sampling.  The  program  uses  pointers  (M(*))  to  indicate  the  S(*) 

Insert  Fig.  1  about  here. 

states  that  can  be  entered  from  *  ,  thus  eliminating  the  nonzero  entries 

in  the  transition  matrix  and  reducing  storage  space. 

Of  special  interest  is  the  method  of  determining  the  state  to  be 

entered  by  a  replication  at  each  transition,  as  in  (4).  MCHAIN  uses  the 

cutpoint  method  developed  in  Fishman  and  Moore  (1981)  whose  execution 

time  is  independent  of  S(J)  ,  the  number  of  states  that  can  be  entered  from  J 

Therefore,  using  this  method  makes  the  cost  of  state  transition  determination 

independent  of  the  number  of  nonzero  entries  in  each  row  of  p  .  MCHAIN 

~n 

uses  the  random  number  generator  GGUBS  in  the  IMSL( 1982)  library,  but 
this  can  be  changed  at  the  user's  discretion. 


MCHAIN  allows  for  three  different  levels  of  output  through  the 


variable  DETAIL.  Setting  DETAIL=2  results  in  a  printout  of  transition 
frequencies  for  each  replication  and  summations  of  these  frequencies 
across  replications.  Setting  DETAIL=1  leads  to  a  printout  only  of  the 
summations  across  replications  and  DETAIL=0  suppresses  all  printed 
output.  In  practice,  we  envision  that  MCHAIN  may  be  used  to  qenerate 
input  for  other  simulation,  as  exemplified  in  Fishman  (1982b).  In  this 
case  a  user  needs  to  alter  MCHAIN  to  store  or  pass  the  frequency  data 
and  may  have  limited  or  no  interest  in  the  printout  of  the  frequency 
data  themselves.  DETAIL'S  options  enable  him  to  specify  his  level 
of  interest. 

MCHAIN  allows  for  two  types  of  simulation  runs,  macrorepl ications 
and  microreplications.  MCHAIN  runs  I  macroreplications  each  consisting 
of  K  microreplications  based  on  parallel  simulation  using  rotation 
sampling.  The  need  for  these  macroreplications  arises  when  estimating 
the  variances  of  estimators.  See  Fishman  (1982b).  In  general,  it  is 
advisable  to  make  K  large  relative  to  I  to  gain  the  benefit  of 
rotation  sampling.  Space  requirements  for  the  arrays  in  MCHAIN  are 
roughly  4  *  (SIZE  +8*NP+2*K+6*  ALL)  . 

Example 

Consider  a  single  server  queueing  model  for  which  the  time  between 
successive  arrivals  are  independent  and  exponentially  distributed  with 
rate  x  ,  service  times  are  independent  and  exponentially  distributed 
with  rate  u>  ,  the  queue  discipline  is  first-come-first-served  and 
there  is  a  finite  capacity  n  .  Corresponding  to  this  formulation  is 
an  embedded  nearest  neighbor  Mar'  "  '•hai  <nose  states  are  the  number  of 


jobs  in  the  system  and  whose  nonzero  transition  probabilities  are 


»01  =  1 


pi,i-l 

II 

i=2,. . . ,n 

pi,i+l 

=  _A_ 

A+w 

i=l ,. . . ,n-l 

(7) 

_  A 

pn,n 

A+u)  * 

We  show  a  run  of  MCHAIN  for  a  first  passage  from  state  a=4  to  state 
b=l 7  for  a  single  macroreplication  containing  15  microreplications  based  on 
rotation  sampling  with  A=.9  ,  w=l  and  n=19  .  The  input  for  MCHAIN 
is  NP=20,  INITIAL=4  ,  ABS0RB=17,  1=1,  K=15,  SEED=1 556203872  ,  SIZE=1000 
and  DETAIL=2  .  The  quantities  ALL  ,  KK  and  N  and  the  K  ,  P  ,  S  and 
SUMS  arrays  are  computed  in  the  sample  driver  program  in  Fig.  2.  Figure  3 

Insert  Figs.  2  and  3  about  here. 

shows  the  output  for  MCHAIN  .  The  entries  in  the  table  are  the  frequencies 
of  transitions  for  all  possible  transition  type  by  microreplication. 
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c  *** 

C  ***  SUBROUTINE  NCHAIN 

C  *** 

C  *** 

C  **♦  THIS  SUBROUTINE  SIMULATES  K  REPLICATIONS  OF  AN  (N+1)-STATE 
C  ♦**  HARKOV  CHAIN  IN  PARALLEL  USING  ROTATION  S  APPLING. 

C  *** 

C  ***  FOR  REFERENCE  SEE: 

C  *** 

C  ***  FISHMAN,  G.  S.  (1982).  "SU PER  EE FIC IENT  SIMULATION  OP  HARKOV 
C  ***  CHAINS  AND  SEMI-MARKOV  PROCESSES,”  REPORT  UNC/ORSA/TR  82/5 
C  ***  AND  "GENERATING  PARALLEL  CORRELATED  TRANSITION  FREQUENCIES 
C  ***  FOR  A  MARKOV  CHAIN,"  REPORT  UNC/ORSA/TR  82/8,  CURRICULUM 
C  IN  OPE  RATIONS  RESEARCH  AND  SYSTEMS  ANALYSIS  ,  UNIVERSITY 

C  **♦  OF  NORTH  CAROLINA  AT  CHAPEL  HILL. 

C  *** 

SUBROUTINE  MCHAIN  (NP ,  IN  ITAL  ,ABSORB  ,P ,Q  ,  H,S  ,SU MS  ,ALL,KK,I,K, 

1  SEED  , SIZE, UU,V, DETAIL, CUTP, ROTATE, CN,KSTAR,KPRIHE, STATE) 

C  *** 

C  ***  DESCRIPTION  OF  VARIABLES: 

C  *** 

C  ***  ABSORB  =  ABSORBING  STATE 

C  ***  ALL  =  SUM  OF  S(J)  FOR  ALL  J 

C  *•*  CN  ( *)  =  HOLDS  TRANSITION  COUNTS  FOR  EACH  TRANSITION  TYPE  AMD 

C  EACH  MIC R ORE  PLICATION;  LAST  COLUMN  CONTAINS  TOTALS 

C  ***  BY  TRANSITION  TYPE  FOR  ALL  MICROREPLICATIONS 

C  ♦♦♦  COLMAX  =  MAXIMUM  NUMBER  CF  COLUMNS  TO  PRINT  PER  OUTPUT  PAGE 

C  ***  CCUNT  =  UNIFORM  DEVIATE  COUNTER 

C  ***  CUT P  (♦)  =  POINTERS  FOR  CUTPOINT  SAMPLING  METHOD 

C  ***  SEE  FISHMAN,  G.  S.  AND  L.  R.  MOORE,  III  (1981). 

C  ***  "SAMPLING  FROM  A  DISCRETE  DISTRIBUTION  WHILE 

C  ♦**  PRESERVING  MONOTCNICITY",  TECHNICAL  FEPORT  81/7, 

C  OPERATIONS  RESEARCH  AND  SYSTEMS  ANALYSIS, 

C  ***  UNIVERSITY  OF  NORTH  CARCLINA  AT  CHAPEL  HILL. 

C  »**  DETAIL  =  DETAIL  FLAG,  VALUES: 

C  ***  0  FOR  NO  PRINTING,  1  FOP  SUMMARY  TOTALS  ONLY, 

C  2  POR  FULL  DETAIL  PRINTING 

C  •**  CSEED  =  RANDOM  NUMBER  GENERATOR  SEED  (DOUBLE  PRECISION) 
c  EP  =  STEP  VARIABLE  FOR  SAMPLING  TRANSITION  TYPE 

C  ***  GGUBS  =  IMSL  SUBROUTINE  FOR  FANDOM  NUMBER  GENERATION 

C  ♦♦♦  I  =  NUMBER  OF  HACROREPLIC ATIC NS 

C  ***  II  =  INDEX  USED  EY  OUTPUT  SECTION 

C  *♦*  IFIRST  =  NUMBER  OF  FIRST  COLUMN  TO  PRINT  CN  CURRENT  OUTPUT  PAGE 
C  ***  II  =  INDEX  USED  EY  OUTPUT  SECTION 

C  IIS  =  S  (J+1) 

C  I  LA  ST  *  NUMBER  OF  LAST  COLUMN  TC  PRINT  ON  CURRENT  OUTPUT  PAGE 

C  ***  IND  =  FLAG,  INITIALLY  SET  TO  0  TO  ALLOW  FIRST  TRANSITION 

C  ***  IF  STARTING  IN  ABSORBING  STATE 

C  ***  INITAL  *  INITIAL  STATE 

C  ***  IP  *  INDEX  USED  BY  OUTPUT  SECTION 


Fig.  1  MCHAIN  FORTRAN  Program 


C  ***  I SEED  =  SAVED  INITIAL  SFED 
C  ***  IX  =  INDEX 

C  ***  J  =  INDEX  FOE  STATES 

C  ***  J1  =  INDEX  FOR  MICROREPLICATIONS 

C  *♦*  JA  =  INDEX  FOR  ALL 

C  ***  JJ  =  INDEX  USED  FOR  INDIRECT  SUBSCRIPTING 

C  ***  K  =  NUMBER  OF  PARALLEL  MICRCREP  LIGATIONS 

C  ***  K1  =  K+1 

C  ***  KA  =  NUMBER  OF  MICROREPLICATIONS  THAT  HAVE 

C  ***  RETURNED  TO  ABSORBING  STATE 

C  **♦  KPPI HE  (*)  =  NUMBER  IN  STATE  AT  END  OF  TRANSITION 

C  **♦  KSTAR  (*)  =  NEXT  KPRIHE  (*) 

C  K1EST  =  TEMPORARY  VARIABLE  USED  IN  OUTPUT  SECTION 

C  *♦*  KTOTAL  =  TOTAL  NUMBER  OF  TRANSITIONS 

C  ***  LCOUNT  =  INDEX  USED  BY  OUTPUT  SECTION  TO  COUNT  OUTPUT  LINES 
C  ***  LI  =  NUMBER  OF  CURRENT  MACRO  REPLICATION 

C  LMAX  =  MAXIMUM  NUMBER  OF  LINES  TO  PRINT  PER  OUTPUT  PAGE 

C  LX  =  CUTPOINT  TO  ACCELERATE  SAMPLING  PROCEDURE 

C  ♦**  M  (*)  =  STATES  THAT  CAN  BE  RFACHED  FROM  J-1 

C  ***  N  =  NUMBER  OF  HIGHEST  STATE  (NP-1) 

C  ♦**  NP  *  NUMBER  OF  STATES 

C  ♦**  NSKIP  =  NUMBER  CF  CUT  PUT  PAGES  TO  BE  PRINTED 

C  ***  CM  =  USED  IN  SETTING  UP  STATE  TRANSITION  PROBABILITIES 

C  ***  P(*)  =  TRANSITION  MATRIX 

C  **♦  Q (*}  =  VECTOR  OF  CUMULATIVE  PROBABILITIES  FOR  EACH  TRANSITION  TYP1 

C  ♦**  R  =  INDEX  FOR  STATES 

C  ***  ROTATE  (♦)  =  ROTATION  INDEX  FOR  SAMPLING  FROM  EACH  STATE 

C  ***  S(J)  =  NUMBER  OF  STATES  THAT  CAN  BE  ENTERED  FROM  J-1 

C  ***  SEED  =  RANDOM  NUMBER  GENERATOR  SEED 

C  ***  SIZE  =  BLOCK  SIZE  FOR  RANDOM  NUMBER  GENERATION 

C  **♦  SJ  =  S(J) 

C  SK  =  SU  MS ( J) 

C  ***  STATE(^)  =  CURRENT  STATE  OCCUPIED  BY  MICROREPLICATION  J1 
C  ***  SUM  S  (J)  =  SUM  OF  S(1)  THRU  S  (J- 1)  ,  AND  SUMS(1)=0 

C  ***  U  =  CURRENT  UNIFORM  DEVIATE 

C  ***  UU(*)  =  CURRENT  UNIFORM  DEVIATE  FOR  TRANSITION  TYPE  ♦ 

C  ***  V(»)  =  UNIFORM  DEVIATE  ARRAY 

C  *** 

INTEGER  CUTP  (ALL)  ,  ABSORB,  ALL,  CN  (KK)  , COUNT,  DETAIL,  KPRIHE  (NP)  ,  I, 

1  1 1 , 1 ND  ,  I N  ITAL,I  SEE  D  ,1  X,  J  ,  J 1 ,  JA  ,  J  J  ,K  ,  KK,  K,  KK  1 ,  K  ,  KK  A,  LI,LX,LCOUi 

2  LMAX, M  (ALL)  , KSTAR  (NT)  ,  NP  ,R  ,S  (NP)  ,  SEED, SIZE,  SJ,  SK, 

3  STATE(K)  ,SUMS(NP) 

INTEGER  CCLMAX, I FIRST, II, I IS, IL AST, IP, ITOTAL, KTEST, KTOTAL, N, NSKIP 
REAL  V  (SIZE) 

DOUBLE  PRECISION  EOT  AT  E  (N  F)  ,  EP,OH,  UU  (NP)  ,  U,  DSEED,P  (ALL)  ,Q  (ALL) 

C  *** 

C  INITIALIZE  DOUBLE  PRECISION  SEED 

C 

CSEEC*SEID 

K1=K+1 
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c  **♦ 

c  **♦  DETERMINE  TRANSITION  PROBABILITIES 

C  ♦** 

DO  100  J=  1  »NP 
SK  =  SUMS  ( J) 

Q(SK+1)  =  P(SM1) 

IF(S(J)  . LT.2)  GO  TO  100 
SJ  =  S(J) 

DO  90  R=2,SJ 

C(SK+R)  *  Q (S K+R-  1)  4  P(SK4R) 

90  CONTINUE 

Q (S  K-f  SJ)  =  1 .0  DO 
100  CONTINUE 
C  *** 

C  ♦**  SET  UP  INDIRECT  INDEXING  POB  TRANSITION  ARRAYS 
C  *** 

DO  200  J=1 ,NP 
SJ=S(J) 

EP=  (1-0/SJ) 

OM=0 

LX  =  SUHS  (J)4l 
DO  200  R= 1 , SJ 

150  CUT P  (SUMS  (  J)  4  R)  =  LX 

IF  (Q  (LX)  .  GT.OH)  GO  TO  195 

LX=LX41 

GO  TO  150 

195  OH-OH+EP 

200  CONTINUE 

DO  950  LI*1,I 

C  *** 

C  *♦*  SET  STARTING  SEED 

C  *♦* 

ISEBD*DSEED 

C  *** 

C  ♦♦♦  SET  FLAG  TO  ALLOW  INITIAL  TRANSITION  FROM  ABSORBING  STATE 
C  *** 

IND=0 

C  *** 

C  GET  ARRAY  OF  RANDOM  DEVIATES 

C  *** 

CALL  GGOES  (IS EED, SIZ E, V) 

C  •** 

C  ***  GET  INITIAL  RANDOM  DEVIATE 

C  *** 

DU(INITAL41)*V(1) 

U*V  (1) 

COUNT-2 

C 

C  *♦*  START  HICBOB  EPL  IC AT  10 NS  IN  INITIAL  STATE 
C  AND  INITIALIZE  TRANSITION  COUNTERS 
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C  *** 

KPBIHE  (INITAL4I )  SK 
DO  210  J1  =  1  #  R 

STATE  (J1)=IN IT AL 
DO  210  JA  s1#  ALL 

CN  ( All* ( J 1-  1)4 JA)  =0 
210  CONTINUE 

DO  215  J=1 #  NP 

FOT  ATE  (J)=0.  0D0 
215  KSTAR  (J)  =0 
KA=0 

C  *♦* 

C  ***  BEGIN  SIHULATION  OF  STATE  TBANSITICNS 
C  •  ** 

220  DO  300  J 1  =  1 ,K 
C  *** 

C  ♦**  DC  NOT  HOVE  FBOH  AESORBING  STATE  EXCEPT  DURING  INITIAL  TRANSITION 
C  *** 

IF  (STATE  (JI).EQ.  ABSORB.  AKD.IND.EC-  1)  GO  TO  300 
I Xs ST  ATE  (J  1)  41 

C  *** 

C  ***  TEST  TO  PREVENT  DIVISION  BX  ZERO 

C  *** 

IF  (K  PRIME  (IX).EQ.O)  GO  TO  300 

C  +♦♦ 

C  ♦**  GET  RANDCH  DEVIATE  FCB  THIS  HICROR EPLICATION 

C  **♦ 

U=UU(IX)  4(80TATE(IX)/KPRlMB(IX)  ) 

IF  (U.GE.  1.0)  U=(J-  1. 

C  *** 

C  ♦**  FIND  CUTPCINT  FOR  SEARCH  OF  SAMPLING  DISTRIBUTION 

C  *** 

LX-  (S  (I X)  *U) 

c  ♦  ** 

C  ♦*♦  DETERMINE  THE  DIB1CTION  OF  TRANSITION 
C  *** 

JJ=CUTP(SUMS  (IX)  4LX4  D 
225  IF(U.LE.Q(JJ)  )  GO  TO  250 

JJ-JJ41 
GO  TO  225 

C  **♦ 

C  *♦*  INCREMENT  THE  AFPEOFBIATE  TRANSITION  COUNTER 

C  •»* 

2  50  CN  <ALL*(J1-1)4JJ)sCN(ALL*(J1-1)4J>?)41 

C  *** 

C  *•*  SIT  THE  STATE  FOB  MICBOR  EPL ICATION  J1 
C  *** 

STATE  (  J 1)  s  R  ( J  J) 

C 

C  ♦**  INCREMENT  THE  TEMPORARY  COUNTER  AND 
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C  •**  THE  ACTIVE  N ICROR  EPLICATION  COUHTER 

C  ♦** 

KSTAR(M(JJ)  -f  1)  *KSTAR(M(JJ)+1)41 
HOT  ATE  (IX)  =ROTATE (IX)41 .0 
300  CONTINUE 

C  *** 

C  ♦**  SET  FLAG  TO  PREVENT  TRANSITIONS  FRCH  ABSORBING  STATE 

C  •** 

IND*  1 

C  ♦** 

C  ***  ADD  NEWLY  ABSORBED  HICROREPLICATIONS  TO  COUNTER 

C  *** 

KA-KA4KSTAR  (ABSORB+  1) 

C  *** 

C  ***  IF  ALL  K  N ICROR EPLICAIIONS  HAVE  BEEN  ABSORBED,  EXIT 

C  *** 

IF(KA.EQ.K)  GC  TO  500 
KSTAR(#ESORB4l) =KA 

C  *** 

C  ♦**  REINITIALIZE  COUNTERS  FOR  NEXT  SET  OF  TRANSITIONS 

C  *♦* 

DO  <100  J=1,NP 

ROT  ATE  (J)  -0.  0D0 
KPRINE  (J)  =KSTA  R  (J) 

KSTAR  (  J)  =0 

IF  (KPRINE  (J)  *  (J-ABSORB-I).EQ.O)  GO  TO  400 

C  *** 

c  ***  GET  NEXT  RANCOR  DEVIATE 
C  *** 

UU(J)  =  V  (COUNT) 

CCU NT=COUNT+ 1 

IF  (COUNT. LE.  SI  ZE)  GO  TO  400 

C  **♦ 

C  •**  IF  NECESSARY,  GET  NEW  ARRAY  OF  RANDOH  DEVIATES 
C  *** 

CALI  GG UBS (DS EEC, SIZE, V) 

COUNT* 1 

400  CONTINUE 

GO  TO  220 
500  SEED=DSEED 

IF  ( DETAIL.  EQ.O)  GC  TC  990 

C  **• 

C  **♦  PRINT  RESULTS  OF  SIEULATION 
C 

W RITE ( 3, 9C00)  LI 

9000  FOR  HAT  ( '  1 '/'  RESULTS  FRCH  HCHAIN  FOR  N A CRO REPLICATION  NO.«,I5 ///) 
550  WRITE(3,90 10)  HP, IN ITAL , ABSORB ,ALL ,K, I, I SEED, SEED, SI ZE 
9010  FORMAT  (///, 

1  '  NO.  CF  STATES  =  ',110//, 

2  •  INITIAL  STATE  *',110//, 
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3 

4 

5 

6 

7 

8 
9 


ABSORBING  STATE 

TOTAL  NO.  OP  (I,J)  PAIRS 

NO.  OF  CORRELATED  HICBORBPLICATIONS 

NO.  OF  INDEPENDENT  HACROREFLICATIOHS 

INITIAL  SEED 

PINAL  SEED 

BLOCKING  FACTOR 


c  * ♦ 1 10//, 
*  *, 110//, 
=*.110//, 
*'.110//, 
=  '.110//, 
=  *,no//> 


c  *** 

C  ***  CALCULATE  TOTALS  ACBOSS  ALL  CHAINS  AND  STORE  IN  CN(*,K1) 
C  *** 

K TOTAL *0 
DC  650  J  A*  1 ,  ALL 
IT01AL  *  0 
DO  600  J1*1 , K 

ITOTAL  *  ITOTAL  4  CN  (ALL*  (J1- 1 )  4JA) 

600  CONTINUE 

CN(ALL*(K1-1)  +  JA)  =  ITOTAL 
K TOTAL  *  KTOIAL  4  CN  (ALL*(  Kl-1)  4JA) 

650  CONTINUE 


C  *** 

C  ***  DETERMINE  OUTPUT  FORMAT  PARAMETERS  (PAGE  WIDTH  AND  LENGTH) 


C  *** 

COLHAX*16 
LMA X=55 
ILAST*  0 

IP  (DETAIL.  NE.1)  GO  TO  680 

C 

C  **♦  PRINT  SUMMARY  TOTALS  ONLY 
C  *** 

WRITE  (3, 2000) 

2000  FORMAT (* 1TRANSITION  *,6X, *RO«*/, 

1  •  FROM  TC*,6X, ‘TOTAL*/, 

2  - - •) 

IP=0 

N=NP-1 
LCOUNTaO 
DO  670  Ja0,N 
IIS-S ( J41 ) 

DO  660  1I=1,IIS 
LC0UNT*LC0UNT41 
IF=IP4  1 

WRITE  (3,2200)  J,  M  (I  P)  ,CN  (ALL*(  Kl-1)  *  IP) 
660  CONTINUE 

IF(LCOUNT.LE.LMAX)  GO  TO  670 
WRITE  (3,2000) 

LCCUNT=0 
670  CONTINUE 
GO  TO  950 

C  •** 


C  **•  DETERMINE  HUMBER  OF  OUTPUT  PAGES 
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c  ♦** 

680  IF(K1.  GT.  COL  MAX)  GO  TO  700 
COLHA  X=K1 
NSK  IP*  1 
GO  TO  750 

700  NSKI P=K1/CCLHAX 

KTEST*  NSKIP*COLN AX 

IP  (  (K  1-KTEST)  .GT.  0)  NSKI  P=  NSKI P*  1 

C  *** 

C  ♦**  PRINT  TRANSITION  COUNTS  IN  TABLE  FORMAT 
C  *** 

75C  DO  900  I 1= 1 , NSKI P 
IFIFST=ILAST+1 
I LAST=I 1 *COLMAX 
IF  (11. EQ- NSK IP)  1LAST=K1 

C  *** 

C  ***  PRINT  PAGE  HEADING 

C  ♦** 

WRITE  (3  ,2  100) 

2100  FORMAT (//* 1TRANSITION  •  ,W2X, *R EPLICATES’ ) 

WRITE  (3,2110)  (IX,IXsIFIBST,ZLAST) 

2110  FORMAT ( *  FROM  TC  *,16(1X,I5)) 

WRITE  (3,2120) 

2120  FORMAT  ('  - *, 

1  i . . —  V) 

IP=0 

N=NE-1 

LCOUNT=0 

C  •  ** 

C  **♦  PRINT  TRANSITION  COUNTS 

C  *** 

DC  850  J=0,N 
1  IS=S  (J-f  1 ) 

DC  800  11=1, IIS 
LC0UMT=L COUNT* 1 
IP=IP*1 

I  A*  AIL*  (IFIBST-  1)  +  IP 
IB=ALL*(ILAST-1)+IP 

WRITE  (3 ,2200)  J,N(IP),  (CN  (IX)  ,  IX=IA,  IB ,  ALL) 

2200  FORMAT  ( *  *  ■ ,  12  ,4  X ,  12  ,  •  -  *  ,1b  (IX  ,15) ) 

800  CONTINUE 

IF  (L COUNT. L E.LM AX)  GO  TO  850 
WRITE  (3,2  100) 

HRITE(3,2  110)  ( IX, IX  =  1  FIR  ST,  ILAST) 

WRITE  (3,2  120) 

ICCU NT*0 

850  CONI1NUE 
900  CONTINUE 

MBIT Z (3,2300)  K1 

2300  FORMAT  (•  0  ***COLUMN  ',15,*  CONTAINS  THE  PCW  TOTALS’) 
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950  WHITB(3,2400)  KTOTAL 

2400  FOR HAT  ( ' 0  •♦♦TOTAL  FOR  ALL  ROUS  IS: 
990  R STORM 
END 
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C  *** 

C  ♦**  THIS  IS  A  SAMPLE  DRIVER  PROGRAM  FOR  SUBROUTINE  HCHAIM 
C  *** 

C  ***  IHIS  DRIVER  PROGRAM  PROVIDES  THE  SETUP  FOR  USE  OF  MCHAIN  WITH 
C  ***  THE  M/M/1/N  QUEUEING  MODEL  WITH  ARRIVAL  RATE  LAN,  SERVICE 
C  ***  RATE  H,  ONE  SERVER  AND  CAPACITY  N.  WITH  THE  EXCEPTION  OF 
C  H,  P,  S  AND  SU  MS,  ALL  OTHER  ARRAYS  ARE  USED  IN  MCHAIN. 

C  *♦* 

C  ***  DESCRIPTION  CF  VARIABLES  - 
C  *** 

C  ***  ABSORB  =  ABSORBING  STATE 

C  ***  ALL  =  SUM  OF  S(J)  FOR  ALL  J 

C  **♦  CN(*)  =  HOLDS  TRANSITION  COUNTS  FOR  EACH  TRANSITION  TYPE  AND 

C  ***  EACH  MICROREPLICATION;  LAST  COLUMN  CONTAINS  TOTALS 

C  ♦  ♦♦  CUTP(*)  =  POINTERS  FOR  CUTPOINT  SAMPLING  METHOD 

C  ***  DETAIL  =  DETAIL  FLAG,  VALUES: 

C  ♦**  0  FOB  NO  PRINTING,  1  FOR  SUMMARY  TOTALS  ONLY, 

C  ***  2  FOR  FOIL  DETAIL  PRINTING 

C  ***  I  =  DESIRED  NUtiBER  OF  INDEPENDENT  HACROR  EPLI  CAT  IONS 

C  ***  INITAL  =  INITIAL  STATE 

C  J  =  INDEX  FOR  STATES 

C  ***  K  =  NUMBER  OF  PARALLEL  MICRCREPLICATIONS/MACROREPLICATION 

C  *♦*  KPRINE  (♦)  =  NUMBER  IN  STATE  AT  END  CF  TRANSITION 

C  ♦  •*  KSTAR  (*)  =  NEXT  KPRINE  (*) 

C  ♦  ♦*  LAM  =  ARRIVAL  RATE  (<W) 

C  ***  M  (*)  =  STATES  THAT  CAN  BE  REACHED  FROM  J- 1 

C  ***  N  =  NUMBER  CF  HIGHEST  STATE  (<50) 

C  •**  NP  =  NUMBER  OF  STATES  f N 4  1  > 

C  ♦**  P{*)  =  TRANSITION  MATRIX 

C  *♦*  Q (*)  *  VECTOR  OF  CUMULATIVE  PROBABILITIES  FOR  EACH  TRANSITION  T 

C  ***  ROTATE!*)  =  ROTATION  INDEX  FOR  SAMPIING  FRCN  EACH  STATE 
C  *♦*  SEED  =  RANDOM  NUMBER  GENERATOR  SEED 

C  ***  SIZE  *  BLOCK  SIZE  FOR  RANDOM  NUMBER  GENERATOR  <<1001) 

C  ***  S<J)  =  NUMBER  OF  STATES  THAT  CAN  BE  ENTERED  FROM  J-1  (<NP) 

C  ♦♦*  SUMS  (J)  »  SUM  OF  S  <1 )  THRU  S(J-1),  AND  SUHS(1)=0  {<100) 

C  ***  DU (*)  =  CURRENT  UNIFORM  DEVIATE  FOR  TRANSITION  TYPE  * 

C  ***  V  (*)  =  UNIFORM  DEVIATE  ARRAY 

C  ***  H  =  SERVICE  RATE 

C  ♦** 

INTEGER  CUTP  (100)  ,  AESOR  B,ALL,CN  (2 500 )  ,DETA IL, INITAL,  I ,  J,  1,  JA, 

1  K,KK, KPRINE  (50)  ,KSTAR(50)  ,  M  (  100)  ,  N,NP, 

2  S  (50)  , SEED, SIZE, STATE  (50)  ,SUNS(50) 

REAL  V  (1000) 

c  m*  00081,8  P8BCISI0N  EOT  AT  E  (50)  ,LAH,  P  ( 100)  ,  Q  (100), UU  (5  0),  H 

C  **•  INITIALIZE  VALUES 

C  *** 

READ  (1,1000)  NP, INITAL, ABSORB, I, K, SEED, SIZE 
1000  FORMAT  (I5/I5/I5/I5/I5/1 12/15) 

WRITE (  3,  70€0)  I  NI TA L .ABSORB 
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7080 

7081 

7082 

1002 

1001 

7085 

C 

c  *•* 
c  *** 


FORMAT!'  INITAL  STATE:  *  ,  15,  *  A  ESORBING  STATE:  *,15) 

WRITE (3,7061)  K,NP 

FOR  HAT  ('  NUMBER  OF  REPLICATES:  ',  IS,  •  NUMBER  OF  STATES:', 15) 

WRITE  (3,7062)  SEED, SIZE 

FORMAT!'  I  NIT  AL  SEED:  '",112,  •  BLOCKING  FACTOR:', 16) 

READ  (1 , 1002)  DETAIL 
FORMAT  (15) 

READ (1,  100 1)  LAN, W 
FORMAT  (2F5.2) 

WRITE  (  3,7085)  LAN, W, DETAIL 

FORMAT  ('  LAM:  ',F5.2,'  W:  ',F5.2,'  DETAIL:  ',13) 

FOR  EACH  STATE,  DETERMINE  THE  NUMBER  OF  STATES  THAT  CAN  BE  ENTERED 
S(1)  =  1 

DO  80  J=2,NP 
S(J)=2 
CONTINUE 

WRITE  (3,2015)  (S  (J)  ,J  =  1  ,NP) 

FORHATC  S(J)  =•  ,2015) 

SUHS(1)  =  0 
ALL-S  ( 1) 

DC  90  J=2,NP 
SUMS  (J)  -ALL 
A1L=AIL4S(J) 

CONTINUE 
KK-ALL*  (K-f  1) 

RESTRICTIONS  ON  PARAMETER  SIZES  BELOW  ARE  FOR  THIS  PROGRAM. 

COMPOTE  TRANSITION  PROBABILITIES 

DC  100  J=2,NP 

P  (SUBS  (J)  41)  *  W/(LAN4W) 

P (SUMS  (J)  42)  *  LAM/(LAM4W) 

CONTINUE 
P(1)  =1  .0D0 

WRITE  (3,2016)  (P(JA) , JA=1, ALL) 

FORMAT  ('  P(*)  *  '  ,  20F5.  2) 

FOR  EACH  STATE  J  DETERMINE  STATES  THAT  CAN  BE  ENTERED 
N  =N  P- 1 

DO  400  0=2, N 

M  (SUMS  (J)  41)  *  J-2 
M(SUHS(  J)  42)  -  J 
CONTINUE 
M(1)  =1 

H  (SUMS  (NP)  41)  =  N-1 

M  (SUHS  (NP)  42)  =  N 

WRITE(  3,  20  17)  (M(JA)  ,JAS1,ALL) 


NUMBER  OF  STATES: ',15) 


2015 


C  *** 

c 

c  **♦ 
c  *** 


2016 

c  *** 
c  *♦* 
c  *♦* 
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f 


2017  POBNAT  |*  «  (*)  =  *,2015) 

C  ♦** 

C  ***  CALI  SUBFOUTINE  M  CHAIN 
C  ♦** 

CALL  NCHAIN (NP,  INITAL,AESOFB,  P,Q,N ,S , SUN S , ALL,K K, I, K, 

1  SEED, SI ZE, UU, V ,DETAl L, CUTP ,BOT ATE, CN, KST AR, KPRIHE, STATE) 

C  *♦* 

C  ♦**  END  OF  PFOGRAH 
C  *** 

STOP 

END 
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estimate  of  0(l/k  )  and  a  computation  time  0(k)  as  k  ®.  This  compares 
favorably  with  the  case  of  independent  replications  wherein  the  estimate 
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example  including  a  sample  driver  program  are  presented  to  illustrate 
how  MCHAIN  works  in  practice. 
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